function res = GetMagnetometerBandwidth(obj, varargin)
    import MeasurementComponent.Signal.*
    
    p=inputParser;
    p.addParameter('nPoint', 50, @(x) x>0);
    p.addParameter('average', 100, @(x) x>0);
    p.addParameter('signal_amplitude', 0.01, @(x) x>0 && x < 0.1);
    p.addParameter('start_frequency', 10.0, @(x) x>0);
    p.addParameter('finish_frequency', 300.0, @(x) x>0);
    p.addParameter('order', 4, @(x) ismember(x, 1:8));
    p.addParameter('bandwidth', 10.19, @(x) x>0);
    p.addParameter('liaIdx', 3, @(x) ismember(x, 1:4));
    p.addParameter('f129', 250.0, @(x) isnumeric(x) && x>0);
    p.addParameter('f131', 75.0, @(x) isnumeric(x) && x>0);
    p.addParameter('fRange', 20, @(x) isnumeric(x) && x>0);
    p.addParameter('isReport', false, @islogical);
    p.addParameter('description', 'normal', @ischar);
    p.parse(varargin{:});

    %%
    if isempty(obj.WorkingPoint)
        error('WorkingPoint cannot be obtained.');
    end
        
    %%
    obj.unLockTransverseField('x'); pause(0.5);
    obj.bx.setDriver( obj.ziSlave.LockInChannel(1).setOscillator(p.Results.liaIdx) );
    obj.bx.setAC('frequency', 10, 'amplitude', p.Results.signal_amplitude, 'offset', obj.WorkingPoint.BiasX.zeropoint);
    liaChannel = obj.ziSlave.LockInChannel(p.Results.liaIdx).setup('order', p.Results.order, 'bandwidth', p.Results.bandwidth, 'phase', 0.0,  'isSinc', true);
    bxSig = Signal('bxAmplitude', liaChannel, 'reference', obj.bx);
    pause(5);
        
    %%
    bxSig.subscribe_signal(ziDemodTime.R_avg, [2, 1, 1]);
    bxSig.subscribe_signal(ziDemodTime.Theta_avg, [2, 1, 2]);
    res.swprData = bxSig.getParameterDependence(obj.ziSlave, ziSweepNode.Osc3Frequency, p.Results.start_frequency, p.Results.finish_frequency, ...
                                          'npoint', p.Results.nPoint, 'averaging_sample', p.Results.average,'isLog', false, 'sincfilter', 1);
    obj.bx.disable();
    liaChannel.setup('isSinc', false);
    obj.LockTransverseField('x');
    
    %%
    res.signal.polarization = 'x';
    res.signal.amplitude = p.Results.signal_amplitude;
    res.signal.start_frequency = p.Results.start_frequency;
    res.signal.finish_frequency = p.Results.finish_frequency;
    res.filter.order = p.Results.order;
    res.filter.bandwidth = p.Results.bandwidth;
    res.filter.sinc = 'on';
    res.timestamp=datetime('now');
    [res.fittingR, res.fittingTheta] = fittingBW(obj, res.swprData, p.Results.signal_amplitude, p.Results.f129, p.Results.f131, p.Results.fRange);

    %%
    obj.MagnetometerBandwidth = res;
    fig = obj.PlotMagnetometerBandwith();
    
    %%
    if p.Results.isReport
        getReporter;
        rpt.NewRecord(['MagnetometerBandwidth_' p.Results.description]);
        rpt.AppendFigure(fig, sprintf('MagnetometerBandwidth %s', p.Results.description));
    end
    %%
    getSession;
    sess.addData('MagnetometerBandwidth', res, p.Results.description);
end

function [fittingR, fittingTheta] = fittingBW(obj, swprData, amp, f129, f131, fRange)
    maxAmp = obj.WorkingPoint.BiasX.slope*amp*obj.WorkingPoint.exportScaleAC/sqrt(2);    
    [fList, dataR] = swprData.getSweeperData(1);
    [xData, yData] = prepareCurveData(fList, dataR);
    
    excludeIndex = find( (f131-fRange<xData & xData<f131+fRange) | (f129-fRange<xData & xData<f129+fRange));
    excludedPoints = excludedata( xData, yData, 'Indices', excludeIndex);
    ft = fittype('a*g/sqrt(x^2+g*g)+c', 'independent', 'x', 'dependent', 'y' );
    opts = fitoptions( 'Method', 'NonlinearLeastSquares' );
    opts.Display = 'Off';
    opts.StartPoint = [maxAmp, min(yData), 20];
    opts.Exclude = excludedPoints;
    
    [fittingR.result, fittingR.gof] = fit( xData, yData, ft, opts );
    fittingR.dcValue = maxAmp;
    fittingR.amplitude = fittingR.result.a;
    fittingR.background = fittingR.result.c;
    fittingR.bandwidth = fittingR.result.g;
    fittingR.exclude = [f129, f131];
    fittingR.excludeRange = fRange;
    fittingR.gain129 = fittingR.result(f129)/amp;
    fittingR.gain131 = fittingR.result(f131)/amp;
    
    %%
    [~, dataTheta] = swprData.getSweeperData(2);
    [xData, yData] = prepareCurveData(fList, dataTheta);
    excludedPoints = excludedata( xData, yData, 'Indices', excludeIndex);
    ft = fittype( 'smoothingspline' );
    opts = fitoptions( 'Method', 'SmoothingSpline' );
    opts.Exclude = excludedPoints;
    [fittingTheta.result, fittingTheta.gof] = fit( xData, yData, ft, opts );
    fittingTheta.exclude = [f129, f131];
    fittingTheta.excludeRange = fRange;
    fittingTheta.phase129 =fittingTheta.result(f129)/pi*180;
    fittingTheta.phase131 =fittingTheta.result(f131)/pi*180;
end
